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Abstract 
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1 Introduction 

We present a simple new approach to the solution of a wide class of spectral and resonance 
problems on infinite domains with regular ends, including those found in the study of quan- 
tum switches, waveguides, and acoustic scatterers. Our algorithm is part analytical and part 
numerical and is essentially a combination of four classical approaches (domain decomposition, 
boundary elements, finite elements and spectral methods) each of which is used in its most 
natural context. 

Our method has three main advantages. Firstly, it allows one to exploit existing library 
software for PDEs to reduce to a minimum the amount of new code which must be written: 
meshing, choice of basis functions, preconditioning and solution of linear systems can all be 
delegated to the library code. Secondly, the method is well adapted to handling resonance 
or scattering problems, where Neumann-to-Dirichlet maps must be continued analytically as 
functions of the spectral parameter. Thirdly, the need to solve boundary value problems for 
many different choices of the spectral parameter is totally removed. 

A brief outline of the solution procedure is as follows. Firstly, the infinite domain on which 
the problem is posed is decomposed (domain decomposition) into a union of a bounded domain 
and a finite number of infinite regular ends, e.g. cylinders. We need to assume that on each 
of these ends, the partial differential equation admits separation of variables. The problem can 
thereby be reduced to a problem on the finite portion of the domain with extra boundary con- 
ditions on the interfaces with the attached ends. This problem can be reduced to a problem 
on the interfaces through the interior 'Neumann-to-Dirichlet' map, whose matrix elements, with 
respect to a specially chosen basis on the interfaces (boundary element method) are calculated in 
terms of traces of 'mixed' Neumann eigenfunctions from the interior domain (spectral method), 
which is probably the main new idea of this article. The Neumann eigenfunctions are them- 
selves computed by a finite element method. The resulting problem on the interface can be 
solved by finding roots of some monotone functions calculated by solving a generalized matrix 
eigenproblem of modest size. This is done by elementary means. 

Apart from the need to be able to separate variables at infinity, the remaining parts of our 
approach work for a wide class of problems with minimal restrictions on the PDE coefficients and 
the boundary conditions. For simplicity of exposition, however, most of this paper is presented 
for the case of the Helmholtz equation in domains with cylindrical ends, or in domains exterior 
to an obstacle. 

Since the 1960s there has been a lot of important theoretical and numerical work done on the 
solution of spectral and scattering problems on the sort of domains we consider. The literature 
is too extensive to review adequately here: we mention just a few trends. 

On the purely mathematical side, attention has recently concentrated on finding conditions 
sufficient for existence of eigenvalues either below the essential spectrum or embedded in it, see 
e.g. [EvLeVal lExSel [DuExl IKrKTj . Additionally, there are some partial results on non-existence 
of eigenvalues in the lower part of the essential spectrum [DaPaj, and on perturbations of 
embedded eigenvalues which turn them into resonances. For an abstract approach to counting 
the eigenvalues, see [Pari IChl IChZw] , 

For numerical treatments, boundary integral and boundary element methods have often 
been employed — see, e.g., Colton and Kress [CoKrj, Brebbia, Telles and Wrobel [BreTeWrj 
for reviews. Others have considered domain truncation methods coupled with the use of 'non- 
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reflecting' boundary conditions — see, e.g., Hagstrom and Keller [HaKe]; Givoli, Patlashenko and 
Keller [GiPaKej . For the particular case of waveguides, special efficient numerical approaches 
can be found in Aslanyan, Parnovski and Vassiliev [AsPaVaj and in Mclver, Linton, Mclver 
and Zhang [ MclLiMclZhj : the scattering matrix approach has also been developed in, e.g., 
Bagraev, Mikhailova, Pavlov, Prokhorovand Yafyasov [BMPPYj and in Semenikhin, Pavlov and 
Ryzhii |SePaRy| . Various related numerical methods and applications to a variety of physically 
motivated special cases and generalisations can be found, e.g., in [Url IBoJol IHaPaGH IPeJol IBoj 
and references therein. 

We do not claim that our approach is more numerically efficient than any of these: what 
is certainly true, however, is that it is much simpler to implement, requiring just a few lines of 
MATLAB code rather than some serious work by expert programmers. 

Some examples of our code are available from the authors' web pages; we encourage the 
reader to download, modify, and use it. 



2 Geometry 

By £1 C d > 2 we denote an unbounded open connected domain with N non-intersecting 
cylindrical ends. Namely, for each 1 < n < N we choose local coordinates (x,y) such that the 
cylinder C° is given by 

C° = {(x,y) : y er°, x > 0}, 

where is the open bounded connected (but not necessarily simply connected) (d — 1)- 
dimensional cross-section of and we suppose that 

TV 

q = q u l_|c°, 

n=l 

where C M d is a connected bounded open set. Further geometrical definitions (interior part 
of f2, cylindrical ends of £1 and interfaces) are introduced in the next section. 

Throughout the paper we assume that the boundary <9S7 is at least piecewise C 1 smooth 
and, for simplicity, satisfies both interior and exterior uniform cone conditions (i.e. dfl does not 
have any cusps). 



3 Spectral problem; boundary conditions; additional notation 

We consider, in $7, the spectral problem 

(3.1) -Au = Xu 

subject to the boundary conditions 



(3.2) Bu :-- 

Here the unit normal n to d£l points outwards 



0. 



We shall refer to Problem (|3TT]), ([3T2]) as "Problem (P)". 

We impose the following restrictions on the boundary coefficients a(-) and &(• 
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(Bi) a and b are piecewise smooth functions on the boundary <9f2; 
(B 2 ) a(-) 2 + 6(-) 2 ^l; 

(B3) on each connected component of the part (I, 00) x dT n of the infinite boundary of each 
cylindrical end C® the functions a and b are constant (but the constants are allowed to 
differ between the cylinders and even between the connected components of the boundary 
of each cylinder). 

To finalize our notation, we shall from now on refer to the sets 

C n :={(x,y) : y x > 1} 

and their union 

N 

C:=[_\C n 

n=l 

as cylindrical ends of 17 and the sets 

r n :={(l,y) <EC n } = {(l,y) : y € T° n } 

and their union 

N 

r := U r n 

n=l 

as interfaces. 

We also define the interior part of O as 

n :=n\c 

and denote by 

r := dn \ r , 

the part of its boundary shared with the boundary of 17, see Figure [I] 

We emphasize once more that according to (B3) the boundary condition (3.2) are constant 
on each connected component of 

N 

dc\r= [J (1,00) x av n . 

71=1 

4 Generalized problem 

As the reader shall easily see later on, all our methods apply, with minimal modifications, to the 
general second order spectral boundary value problem 

(4.1) — div(p(-)grad u) + q()u = Xu 

(4.2) (a(-)« + 6(-)n-p(-)gradu(.))| en = 0. 
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Figure 1: Example of a domain with cylindrical ends. 



with p being a sufficiently smooth positive definite [d x <i)-matrix valued function, and q being 
an Lcxj(R ) scalar potential. Then, in addition to assumptions (Bi)-(B3), we also require that 
p(-) and q(-) be constant on each cylindrical end C° (but the constants are allowed to differ 
between the cylinders). 

Although one can easily formulate all the results and algorithms for the general case (4.1), 



(4.2), we deliberately describe only Problem (P) for the Laplacian in order to preserve the 
simplicity and clarity of exposition. 



5 Rigorous problem statement 

We always understand the spectral problem (P) in the variational sense. Namely, let 



{x G <9S7 : b(x) = 0} be the part of the boundary on which (3.2) is the Dirichlet condition. 
Let Cg° (Q) := {v E C°°(0) : supp v n d D VL = 0} be the set of all infinitely differentiable 
functions on supported away from dp^l- We define the sesquilinear form 

(5.1) t[u,v}:= jvu-Wv+ j °^uv 

n an\d D n 

first on Cq°(0), and then extend it to the closure Hl(tt) of Cq°(Q) with respect to the norm 



associated with (5.1 ) 



From now on, we shall understand Problem (P) as the spectral problem t[u,v] = X(u,v), 
where (•,•) denotes the scalar product in Z^H), with the form domain in Hl(Q). In particular 
A is an eigenvalue of Problem (P) if there exists a nontrivial function u G Hl(Q) such that 
t(u,v) = X(u,v) for all v G Hi (SI). 
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6 Transversal problem and essential spectrum 



We need to consider the following transversal problem on the joint interface T: 

dw s 



(6.1) 



KW 



a{-)w + b{-) 



dv 



0: 



ar 



here A^_i is the (d — l)-dimensional Laplacian and is the derivative with respect to the 
external normal v to dT (which of course coincides here with the normal n to <9S7). We recall 
that a and b are constant on each connected component of dT. 

This problem has a rigorous formulation in terms of sesquilinear forms analogous to that 



introduced in (5.1) but with 0, replaced by T throughout. The assumptions (Bi), (B2) and 



(B3) ensure that the problem (6.1) has a purely discrete spectrum of eigenvalues («»• 



00 

■J J j=l' 

" ' 00 

3=1 



enumerated in increasing order and repeated according to multiplicity; we denote by (wj(y))j 
the corresponding eigenfunctions normalized in ^(T). We also set kq := —00. 

It is well known that the essential spectrum of Problem (P) is then given by [ki,+oo) and, 
moreover, that its multiplicity at a point A > n\ is equal to the number of Kj which are less 
than or equal to A. We shall often refer to the number Kj as the j-th threshold. 

For each J > 1 we define 

Pj := Span{^ : j > J} , 

set Po := ^(T), and denote by Vj the projection in i^CO onto Pj (J > 0). We also define 

Qj:=I- Vj. 



7 Eigenvalue problem — reduction to the interface 

We define formally two Neumann to Dirichlet operators, for the interior part S7o ar| d cylindrical 
ends C. 

Firstly, the interior Neumann to Dirichlet operator 1Z\ := g 1— > v\ r is the map which maps 
any suitable function g defined on T to the trace on T of the solution of the boundary value 
problem 

dv 

(7.1) - Av = Xv in n , — = g on T, Bv = on T 

on 

(where by n we always mean the unit normal to T directed outwards from Oo) provided A is not 
an eigenvalue for the corresponding Neumann problem 

dU 

(7.2) -AU = XU in O , — — = on T, BU = on T . 

on 

Secondly, the exterior Neumann to Dirichlet operator TZ^ := g 1— > V\ T is the map which 
maps any suitable function g defined on T to the trace on T of the solution of the exterior 
boundary value problem 

dV 

(7.3) - AV = XV in C, -— = g on T, BV = on dC\T 

on 

(the normal n on T, as before, points outwards from f2o along the positive x-axis of each 
cylindrical end, and therefore towards C). 
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Suppose that u is an eigenfunction of Problem (P) with eigenvalue A. It is well known that 
u is infinitely differentiable at any interior point of £1. In particular, therefore, u and its gradient 
restrict to the interface V as smooth functions, and so we may define 



(7.4) 

Thus / and g satisfy 



du 
dn 



f = -K c x g = TZ° x g, 



so that a necessary condition for A to be an eigenvalue of Problem (P) is that a 
eigenvalue of the operator pencil problem 



-1 be an 



(7.5) 



(aTZ c x -TZ° x )g = 0. 



(A number ctq is called an eigenvalue of a general operator pencil aA + B if zero is an eigenvalue 
of the operator cfqA + B.) 



8 Problem on a cylindrical end — the Neumann to Dirichlet op- 
erator 



Suppose that V is a solution of the exterior problem (7.3) for some A (not necessarily an 
eigenvalue of Problem (P)). A formal expansion of V in terms of the eigenfunctions Wj of the 
transversal problem ( |6.1[ ), 

oo 

V(x,y) = ^2c j (x)w j (y) 

3=1 

with account of (17.31) yields for the functions Cj the ordinary differential boundary value problems 



(8.1) 



d 2 Cj(x) 
dx 2 



+ (A — Kj)cj(x) =0 for x > 1 



dcj(x) 



dx 



x=l 



where (•, •) denotes the scalar product in ^(r). 
The two solutions of (8.1) are 

(8.2) 



cf(x) 



T^^=exp(±(x-l)v^A) . 

yj Kj — A V / 



If we additionally require that the solution V G H;(C) and assume that, for some in- 
teger J > 0, A G [kj, Kj+i), then the boundary Neumann data g should belong to the 

— 1/2 

space QjH* (T) (in particular, we need (g,Wj) = for j < J) and choose only the so- 
lutions Cj(x) := cj(x) in order to to exclude terms growing or oscillating at x = +oo. Then 

V(l,y) = Y^j=j+i c j(^) w j(y) ^ H* Cn an d we have therefore defined the non-negative 
self-adjoint operators 



K L X : QjH, 1/2 (r) -» QjH: /z (T) for A G [kj, k j+ i) 



1/2, 
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acting as 
(8.3) 



^9= £ 



j=J+i 



(g,Wj)wj. 



It is clear from (8.3) and from the construction that the quadratic form (JZ c x g,g)r is mono- 
tone increasing in A on each interval [kj,kj + i). 



9 Problem on an interior domain — the Neumann to Dirichlet 
operator 

Let g\ and g2 be functions defined on T and suppose that v\ and v 2 are the corresponding 



solutions of the boundary value problem (7.1), for A which is not an eigenvalue of (7.2). A 
formal integration by parts shows that 



(9.1) 



(K9i,92) = (gi,K° x g 2 ) = (Vv l ,Vv 2 ) - \{v 1 ,v 2 ) 



This allows us to define a quadratic form 



t x [g) := (K° x g,g) = (Vv,Vv) - X(v,v) 



(where v solves (7.1)) associated with 1Z\ and hence to regard 1Z\ either as a self-adjoint 
operator on a domain in L 2 (T) or, under suitable smoothness restrictions on the boundary dflo, 
as a map on a scale of Sobolev spaces from H s - l / 2 (T) to H S+1 / 2 (T) (see (Ag) [PaB I5i]). 
We emphasize that we impose a condition that A is not an eigenvalue of ( |7.2[ ) — for ways of 
bypassing this condition see [Si]. 

For a fixed g E L 2 (T), the form x\[g] is known to be monotone increasing with respect to A 



in any interval not containing the eigenvalues of (7.2), as the following simple argument (found 



e.g. in [Fr]) shows. Let g be fixed and let v' denote the derivative with respect to A of the 



solution v of (7.1). Then v' solves the boundary value problem 



(9.2) 



Av' = Xv' + v in Q , dv'/dn\ r = , Bv'\ 



0. 



Multiplying (9.2) by v 1 and integrating by parts we obtain 

1 d 



(9.3) 



2 dX 



tx[g) 



d 
dX 



! dv 
' dv 



{v,v)n > 0, 



which proves the claim. Exactly the same argument can be used to prove the monotonicity in 
A of the form of the operator 1Z X , which we established at the end of the previous section by 



direct inspection of (8.3) 



10 Problem on an interior domain — basis expansions 



In § |8j we constructed an explicit representation (8.3) for the Neumann to Dirichlet operator 
1Z C X associated with the problem on a cylindrical end. The aim of this section is to do the same 



Eigenvalues for domains with infinite regular ends 



M Levitin and M Marietta 



August 2008 (final) 



Page 10 



for the interior Neumann to Dirichlet map H\, and we want to emphasize that the following 
trick is in fact the focal point of this paper. 

In the setting of problem (7.1), let us choose an arbitrary basis {(ftk}™^ in I^CH (which 
is the domain of the quadratic form r^)- Let us also denote by the solution v of (7.1) with 
g = 4> k , i.e. of 



(10.1) 



\®k in SI 



o, 



dn 



b k on T, B<S> k = on r . 



We also denote by /x m the eigenvalues and by U m the corresponding eigenfunctions of the 



homogeneous problem (7.2), i.e 
(10.2) 



— AU m — M m U m in 



dn 



on T, BU m = on r , 



Our aim is to find explicit expressions for the matrix elements R 



By (9.1 ) we have 
(10.3) 



We now use the fact that the eigenfunctions U m form a basis in the Hilbert space HI(£Iq) 
We can thus expand &k in this basis as 

oo 
m=l 



Thus , expanding the right-hand side of (10.3) with the use of two elementary integrations by 
parts, we obtain 



(10.4) 



m=l 

oo 



j^0"m - X){&k,U m ){U m ,$t) = Y — TT^kifmlrJr • (0m|r,&)r- 

m=l ^ m A 



m=l 



11 Embedded eigenvalues and orthogonality conditions 

We have already established formally in § [7] that if A is an eigenvalue of Problem (P) then 
a = — 1 is an eigenvalue of the pencil ( |7.5[ ) on the interface T. In view of rigorous definitions of 
the exterior Neumann to Dirichlet operator, we want now to consider, for A G [kj, Kj+i), the 
pencil 

(11.1) Ax(a):=aK c x -QjK° x . 

acting in a weak sense in Qjl^fi)- 

For each fixed A (which plays the role of a parameter rather than spectral parameter) A\(cr) 
has (T-spectrum spec(^.A), and a necessary condition for A to be an eigenvalue of Problem (P) 
is that —1 G spec(^4.A)- The corresponding eigenfunction g 6 QjZ^CH is the normal derivative 
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on T of the eigenfunction u of Problem (P) corresponding to A. Also, if A £ [kj, Kj+i) for some 
J > then we have seen in § [8] that g must additionally satisfy the orthogonality condition 
VjK\g = (J — Qj)U° x g = 0; in other words TZ° x g G QjL 2 {n). 

Since QjL2^l) is an invariant space of the operator 1Z X , we propose the following strategy 
for finding embedded eigenvalues: 

• first find those "suspicious points" A* for which a = — 1 is an eigenvalue of the operator 
pencil flll.ip on the interface T; 

• if these "suspicious" A* lie in [kj,kj + \) for some J > 0, check the orthogonality condition 
VjV\ g = for the corresponding pencil eigenfunction g G QjL2(T). 

Assume that we wish to find the eigenvalues in the interval [kj, «j+i). We start by reducing 
the spectral problem for the pencil ( | 11 . 1[ ) on the space QjLzfi) to a spectral problem for an 
infinite matrix pencil. In § [8] we developed a representation of the operator 1Z C X in terms of 
the transverse eigenfunctions wf, in §[Io]we developed a representation for TZ° X in terms of an 
arbitrary basis {(pk}kLi in £2 CO- We now choose as our basis of QjL 2 {T) the set {4>k}kLi 
given by cf) k = w k+ j, k G N. Consequently we obtain from (10.4) that 

(11.2) R ki := (TZ x w k+ j,Wi + j)r = V] r(wk+j,U m \ r )r ■ (U m \r,wi +J ) r . 

— 1 fJ>m ~ * 
m=l 

We simplify the expression for the infinite matrix R(A) = {Rke)i<k,e<oo by introducing some 
notation. Define an infinite matrix S with entries {Sk m )i <k m<00 defined by 

(11-3) S km = {wk, U m \r) r , 

an infinite diagonal matrix D(A) with diagonal entries 

(11.4) D kk (X) = (fA k - A)" 1 , l<A;<oo, 
and an infinite diagonal matrix T(A) with diagonal entries 

(11.5) T kk (\) = ^=L=, l<A:<oo. 

For a general matrix M we will use the notation M a: f, jC: d to denote the submatrix 

(11.6) (Mij) a <i< b , 

c<j<d 

where b = 00 and d = 00 are allowed. 



The matrix R whose entries are defined in (11.2) can now be written as 
(11.7) R(A) = S 

J+l:oo,l:oo 

D(A) (Sj+i :00 ,i :00 )* . 

In a similar way, we can define the matrix elements associated with 1Z X with respect to the 
same basis using (8.3): 1Z C X has a diagonal matrix representation T(A) in this basis. We thus 
obtain the following reduction of the problem from QjL2(T) to £ 2 (N). 
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Theorem 11.1. A number A G (which is assumed not to be an eigenvalue of the 

auxilliary Neumann problem (10.2)) is an eigenvalue of Problem (P) if and only if —1 is a 
a-eigenvalue of the matrix pencil 



(11.8) 



crR(A) - T(A) = crSj + i :0Oi i :0O D(A) (Sj+i :00i i :00 )* - T(A) 



and, if J > 0, the correspoding eigenvector c of (11.8) satisfies the orthogonality condition 

(11.9) Si :Jil:00 D(A) (S J+ i :oo ,i :oo )* c = 0. 



Proof. Let A G be an eigenvalue of Problem (P), let u(x,y) G be a corre- 







sponding eigenfunction, and let its traces / and g be defined by (7.4). Then the restriction of u 
to the cylindrical ends C is a square integrable solution of (3.1) and therefore Vjf = Vjg 
(or Qj9 = g) and — 7£?p = —T^Qjg = f, with 7£Ca given by (8.3). Similarly, from the 
restriction on the interior part Oq, we obtain that / = 7£°g which immediately leads to 



(11.10) 



A x (-l)g =(• 



-ft? 



Qjft°)3 = 



and additionally implies the orthogonality condition 
(11.11) VjK° x g = 0. 



If we now write g = Y^k=i c k w k+J with c^ = (g,Wk+j)v and use the definitions (11.2)- 
11.7), then by elementary matrix manipulations (11.10) turns into (— R(A) — T(A))c = and 
TTTTlh becomes rfTOb. □ 



12 More on eigenvalues of the operator pencil (11.1 ) 



Before we proceed to descrtibing in detail the numerical procedure which implements Theo- 
rem 1 1 1 . 1 1 we want to establish some more facts about the eigenvalues <Jj(X) of the operator 
pencil A\(a) = cfTZ c x — Q,jH\ (for A G [kj, kj+i)). We enumerate the pencil eigenvalues <x,(A) 
in decreasing order: o"i(A) > ^(A) > 

Lemma 12.1. The a -eigenvalues of A\(a) are monotone increasing functions of X on each 
interval [Ai,A2] C [kj,kj+i) which does not contain any eigenvalues (j, m of the homogeneous 



Neumann problem (10.2) 



Proof. We established in §§ [8] and [9] that each of the operators 1Z C X and QjT^I is monotone 
increasing (in the sense of corrresponding forms) on each interval satisfying the conditions of 
the Lemma. Thus, for each fixed negative value of a, the quadratic form of the pencil A\(a) 
is monotone decresing in A, and, as 1Z C X > 0, by the inverse function theorem the eigenvalues 



cr(A) increase as functions of A. 



□ 



We also want to estimate how many pencil eigenvalue curves <j(A) can cross the line a = — 1 
as A increases in a given interval (or, equivalently, to estimate a priori how many eigenvalues of 
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Problem (P) can be located in a given interval). We need to introduce some additional notation. 
Let v m denote the eigenvalues of the interior homogeneous Dirichlet problem: 



(12.1) 



AU = uU in a 



o- 



U = on T, BU = on T . 



(This problem is similar to (10.2) with the only difference being that an auxilliary Neumann 
condition on T is replaced by the Dirichlet one.). 

Lemma 12.2. Let A2 G \ {fJ> m }, snd let Ai > max{ftj, max{// m : [i m < A2}}. 

(In other words, we choose an interval [Ai,A2] in such a manner that it does not contain any 
thresholds or any Neumann eigenvalues jjL m .) Then the number of eigenvalues of Problem (P) 
in the interval [Ai, A2] does not exceed 



(12.2) 



K := #{m : [i m < A 2 } - #{m : v m < A 2 } . 



Proof. By Lemma 12.1 the curves Cj(X) are monotone in A G [Ai,A2]. We wish to estimate 
from above the total possible number of times these curves intersect the straight line o = — 1. 
Obviously, this does not exceed the number of curves which are above —1 at the right-hand side 
of the interval [Ai, A 2 ], namely K' = #{j : -1 < aj(A 2 )}. We shall now show that K' < K, 



where K is given by (12.2) 



Note that any A* for which a particular curve Cj(X) crosses a 
ira 
Thus, 

K' < #{positive eigenvalues of - TZ c k2 - QjTl\ 2 } 
= #{negative eigenvalues of TZ C A + QjT^a,} . 



1 is at the same time 



characterised by the fact that zero is in the spectrum of the operator A\+(—l) = —IZ^—QjIZ 



(12.3) 



Since is non-negative for any A, we can majorize the right-hand side of (12.3) by omitting 
it there (and also dropping the projector Qj). Thus, 



(12.4) 



K' < #{negative eigenvalues of 1Z\ } 



But by [Sa] (which is in turn a generalisation of [Frj, see also |Ag|) there is a known relation 
between the number of negative eigenvalues of an interior Neumann-to-Dirichlet operator, and 
the differences of the Neumann and Dirichlet counting functions: 



(12.5) ^{negative eigenvalues of Tl\ 2 } = #{m : p, m < A 2 } - #{m : v m < A 2 } = K . 



□ 



Remark 12.3. The reader may notice that the estimate (12.2) does not in fact use Ai and is 



therefore just an estimate on the counting function of the eigenvalues of Problem (P). One can 
easily improve it if one knowns an estimate from below on the number of negative eigenvalues 
of the operator 1Z c Ai + QjTZ° Ai , say 

#{ negative eigenvalues of lZ Al + QjT^Ail ^ ^1 



(e.g. obtained numerically). Then in the statement of Lemma 12.2 one can replace K by 
K — K\. Further improvements are also possible by introducing the counting functions of a 
generalised Neumann problem instead of the ordinary Neumann problem, but we do not discuss 
them here. 
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13 Numerical algorithm 



Theorem 11.1 gives the characterisation of the eigenvalues of Problem (P) which forms the 
basis of our simple numerical approach to their calculation. We emphasize that this approach 
depends only on having software capable of the following: 

• finding eigenvalues and eigenfunctions for problems on bounded domains with mixed 
boundary conditions (e.g. MATLAB PDE Toolbox, FEMLAB); 

• performing quadratures on boundaries and cross-sections of domains; 

• finding the first few eigenvalues of a matrix pencil; 

• finding the zeros of a monotone real-valued function of a real variable. 
Our algorithm in its basic form is as follows: 

Main Algorithm. 

We remind the reader that the \i m are the eigenvalues of our PDE on £Iq with Neumann 
boundary conditions on the interface T (problem (7.2); see also Fig. [Ip; the U m are the 



corresponding eigenfunctions; and the Wk are the eigenfunctions of the interface problem (6.1) 
with corresponding eigenvalues k^. In two dimensions, the interface problem ( |6.1[ ) becomes an 
ODE problem and the Wk will be explicit trigonometric or hyperbolic functions. 

1. Choose a sufficiently large JV6N and calculate numerical approximations to the Neumann 
eigenvalues (/z m )^ =1 and their eigenfunctions (t7 m )^ =1 for the inner Neumann problem 



(7.2), using an appropriate method - e.g. a standard finite element package. 

2. Choose M G N such that the traces of the eigenfunctions (£/ m )m=i on the interface V are 
well approximated by linear combinations of the first M eigenfunctions (v)j)jLi of the 



interface problem (6.1). Since the convergence of such expansions is usually exponential 



it is rare to need M > 20. 



3. Calculate the first M eigenvalues (Kj)^ =l and eigenfunctions (•Wj)jL 1 of the interface problem 



(6.1) — either analytically in terms of sines and cosines, as in the case of the two- 
dimensional examples treated here, or numerically in the case of general higher-dimensional 
examples. 



Calculate the entries 

r 



(wk, U m \r) r = j Wk U m \ T : 



fork = 1, ... , M, m = 1, . . . , N of the matrix S (see ( 11. 3[ ) ): these are the inner products 



of the interface eigenfunctions Wk with the traces on the interface of the eigenfunctions 



U m of problem (7.2) 



5. Fix J, < J <C M - we shall look for eigenvalues A € [kj, kj+i). Recall that we are using 
the convention kq = — oo, so for eigenvalues below the essential spectrum we choose 
J = 0. No calculations are carried out at this step. 
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6. Calculate the first N Dirichlet eigenvalues of the problem ( 12.1 ) and use the estimates in § [J2| 

to obtain an upper bound K on the number of eigenvalues of Problem (P) in [kj, kj+i). 
Note: K < M. 

7. Define the functions giving the elements D^k, k = 1, . . . , N, of the diagonal matrix D(A) 

and the elements Tkk, k = 1, . . . , M, of the diagonal matrix T(A), by 



■Djfefc(A) = {fjbk ~ A) 1 ; T fcfc (A) = (n k - A) 



-1/2 



Here the kj, are the eigenvalues belonging to the interface eigenfunctions Wf., calculated 
in Step 3. No calculations are carried out at this step. 

8. Define a procedure which takes a value A G [kj,kj+i) and returns the largest K o- 

eigenvalues 

^i(A) ><r 2 (A) > ■■■>a K (X) 
(out of the total of M — J cr-eigenvalues) of the matrix pencil 

(13.1) ffSj+l:M,l:NDl:N,l:NW (Sj+1:M,1-.n) ~ T J+1:M,J+1:m(^) ■ 

K is given in Step 6. No calculations are carried out at this step. 

9. Using the procedure of Step 8 and an appropriate numerical rootfinding algorithm, solve (if 

possible) the K nonlinear equations 

(13.2) o-j(A) = -l, j = l,...,K. 

Note that care has to be taken if any of the interior Neumann eigenvalues p, m lies in 
[k j, kj + i). Away from the \i m , the crj are monotone functions of X and so the rootfinding 
is, in principle, easy. 



Denote by Xjj the solution of (13.2) and by c^j the eigenvector of the pencil (13.1) 
with X = Xj^j, a = — 1. 

10. If J > 0, check the orthogonality conditions 

(13.3) Si:J i i : jvDi : jv,1:Jv(A) {Sj + i : m,1:n) Cj,J = 0. 

(Of course, these can never be checked exactly: at most one can check that they are 
'almost' satisfied and that there is therefore either an embedded eigenvalue or a point of 
spectral concentration near Xj t j.) 

□ 

In the absence of error from the numerical approximation of the eigenvalues and eigenfunc- 
tions, the accuracy of this algorithm is determined by the choice of N at the very first step. 
(The choice of M in step 2 is generally less problematic as the transversal eigenfunctions provide 
rapidly converging approximations to smooth functions on the interface.) N is the number of 
terms at which the sum in the expression (10.4) for Rki(X) is truncated, and unfortunately 
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this series can be rather slowly convergent. We accelerate the convergence with a cheap trick. 
Suppose that for some fixed Ao we already know R^Xq). Then (10.4) yields 

A-A 



(13.4) R M {\) - R M (X ) = 



\ {fJ-m — A)(^ m — Aq) 



m r r 



m r, 



The series in (13.4) converges much more rapidly than the series in (10.4). Moreover the 
numbers Rke(Xo) are easily calculated from ( 10. 3[ ) (with A = Ao) by solving the boundary value 
problems ( 10. 1[ ) numerically (for A = Ao) using any reasonable discretization scheme for elliptic 
PDEs. Our Main Algorithm is then modified by changing steps 8 and 10 as follows: 



8'. Let R° denote the matrix with entries Rke{Xo) defined in (11.2) (with A replaced by Aq). 



Then the matrix pencil (j 1 3 . 1|) is replaced by 
(13.5) 



a 



R 



J+l:M,J+l:M + $J+l:M,l:N {^>1:N,1:N (A) - Di^i^Ao)} (S J+1:M,V.n)* 

-Tj + l : M,J+l:M(A). 



10'. The orthogonality condition (13.3) is replaced by 

(13.6) [Rl-J;0,l:7V + Sl:J,l:JV {Dl:JV,l:Jv(^) ~ Dl:iV,l:7V (^o)} (Sj+l:M,l:7v)*] Cj,J = 0. 

Note that R^j.q i-at ma kes sense: in view of our choice of basis functions = Wk+j, 
its (k,j) element is the numerical approximation to (lZ\ Q Wk, Wj + j)r- 

The convergence acceleration trick can be repeated if one is prepared to evaluate the matrix 
elements R^e numerically for several different A. However for the numerical examples which we 
examined here, one application was always sufficient. 



13.1 Example 1: Bent waveguide 

We consider the Laplacian in a waveguide in M 2 of width 1, bent through an angle ir/A as 
indicated in Fig. [2] (the inner radius of curvature of the bend is 1, the outer radius 2). On the 
side of the waveguide with the smaller radius of curvature on the bend, Dirichlet conditions are 
imposed; on the other side the boundary conditions are of Neumann type. 

We know that the essential spectrum is determined by the spectrum of the problem on the 
cylindrical ends only, and is therefore the same as for the unbent waveguide: 

Cess = [71" 2 /4, +Oo). 

It may also be shown that for the unbent waveguide the spectrum is purely absolutely continuous. 
However this changes when the waveguide is bent. Krejcink and Krfz [KrKrj show that bending 
the waveguide in the direction of the Dirichlet boundary condition will cause an eigenvalue to 



appear below the essential spectrum. Using the method described in § 13 this eigenvalue can 
be located quite precisely. For this purpose the inner domain S7o is taken to be the part of the 
domain between the two dashed lines in Fig. [2] The basis on the cut-off boundary T of S7o 
(the part of dflo composed of the two dashed lines in Fig. [2]) is defined as follows. Let s be 
the local transverse coordinate on T measured with s = on the 'Dirichlet' side. We denote by 
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width = 1 




Figure 2: Bent Waveguide 



02j+i, j = 0, 1,2, . . ., the basis functions supported on the vertical portion of T, and by (f>2j, 
j = 1,2,..., the basis supported on the part of T which is at angle 7r/4 to the horizontal. In 
particular, 4>2j is zero on the vertical part of T and 4>2j+i is zero on the part of T which is at 
angle tt/4 to the horizontal. On their supports, the tfij are given by 



(t>j{s) 



sin(jvrs/2) (j odd), 

sin((j — l)7rs/2) {j even). 



Table [I] shows the results of computations at four different levels of accuracy. 



accuracy 



eigenvalue found 



silly one mesh refinement 

low three mesh refinements 

medium four mesh refinements 

high five mesh refinements 



sum over Xj < 10 in ( 10.4 1 
sum over Xj < 50 in ( 10.4 1 



sum over < 100 in ( 1(14 
sum over Xj < 200 in (IU~4 



none found 
2.3461 
2.3459 
2.3454 



Table 1: Levels of accuracy and eigenvalue found below the essential spectrum, caused by 
bending the waveguide. 



14 Generalization to resonances 



The algorithm in § 13 can be generalized to the calculation of resonances. There are two steps 
in this generalization. 

Firstly, resonances are generally complex. (When a resonance becomes real, it usually be- 
comes either an imbedded eigenvalue, which we have already discussed, or a point of spectral 
concentration — a local maximum of the derivative of one of the spectral measures associated 
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with the operator.) Since resonances are complex, the condition A G [kj, Kj+i) introduced in 
§[8] for the calculation of imbedded eigenvalues no longer makes sense. Of course the equation 
— AV = XV can still be solved in the cylindrical ends, following the procedure i n § [8j and the 
solution expressed in the form 



V(x,y) 



oo 

i=i 



c j (x)w j (y). 



The functions Cj are now chosen to be 



(9,Wj) 



exp 



((x-i)v^a) . 



(cf. (8.2)). The reason for such a choice is that for calculation of resonances we drop the 
condition V £ Hl(C) and we explicitly choose the branch of the square root such that Cj 
£2(0,00). Thus, if we seek resonances (for definiteness, in the lower half-plane ImA < 0), we 
choose the branch of the square root such that Re(y / Kj — A) > 0. 



The expression (8.3) remains formally unchanged: however one must choose the branch of 



the square root as above. A similar remark applies to the matrix T defined in (11.5). With this 



proviso the matrix pencil (13.5) is formally unchanged, though J must be assigned the value 0. 

With these modifications, the resonance approximations are the values of A such that some 
cr-eigenvalue crj(A) of the pencil (13.5) satisfies (Tj(A) = 1. However calculating the resonances 
in this way is impractical. Even if the aj are analytic functions of A, which is not guaranteed, 
it is difficult to order the numerically calculated values of the aj to correspond to an analytic 
ordering. 

The simple approach which we used here was to calculate the condition number of the pencil 

(14.1) [R?:Af,l:M + $1:M,1:N {Dl : jV,l:Ar(A) — Di : jv,l:Jv(Ao)} (Sl:Af,l:Jv)*] _ Ti : m,1 : m(A). 

as a function of A and make a contour plot. Local maxima are 'suspicious points' as the 
condition number will be infinite at a resonance. (The converse, however, is not true: not all 
local maxima are resonances.) Once suspicious points are located approximately one can 'zoom 
in' and examine them in more detail. 



14.1 Example 2: obstructed waveguide 

This problem is considered in detail by Aslanyan, Parnovski and Vassiliev [AsPaVa]: a waveguide 
in R 2 in the strip |y| < 1, obstructed by a plane obstacle O — see Fig. [3] For the particular 
experiments here, we consider an obstacle which is symmetric about the y-axis. The domain 
can then be reduced to the strip 

{(x,y) e R 2 : x > 0, \y\ < 1}\0. 

Following Aslanyan et al. we take as the obstacle a disc of centre (0,5) and radius R, with 

\S\+R< 1: 

= {(x,y) : x 2 + (y - 5) 2 < R 2 }. 
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1 


c 












-1 



X 



Figure 3: Waveguide obstructed by an obstacle 

When 5 = the problem has an eigenvalue imbedded in the essential spectrum. When 5 is 
perturbed from zero, the eigenvalue evolves into a resonance, analytically as a function of S (see 
[AsPa Va]). For the numerical experiments we truncated the strip at x = 1: thus, we took as 
our domain Qo the set 

{(x, y) G M 2 : < x < 1, \y\<l}\0. 

Table [2] shows the results for various values of the numerical experiment parameters. Each is 
quoted for four 'levels of accuracy' defined in Table [T] 



5 




i?=0.3 




e 




R = 0.5 




eigenvalue 


accuracy 


genvalue 


accuracy 




1.50866 




S 


1 


39806 




S 




1.50499 




L 


1 


39157 




L 


0.0 


1.50489 




M 


1 


39139 




M 




1.50486 




H 


1 


39134 




M 




1.5048 




A 


1 


3913 




A 




1.5120 4 


10- 4 i 


S 


1 


4043 4 


8 x 10 4 i 


S 




1.5080 4 


10~ 4 i 


L 


1 


3981 4 


9 x 10 4 i 


L 


0.1 


1.5079 4 


10~ 4 i 


M 


1 


3979 4 


9 x 10 4 i 


M 




1.5078 4 


10~ 4 i 


H 


1 


3979 4 


9 x 10 4 i 


H 




1.5102 4 


10~ 4 i 


A 


1 


3998 4 


9 x 10 4 i 


A 




1.5203 4 


4 x 10~ 4 i 


S 


1 


4236 4 


3.46 x 10 3 i 


S 




1.5167 4 


5 x 10 4 i 


L 


1 


4180 4 


3.89 x 10 3 i 


L 


0.2 


1.5165 4 


5 x 10 4 i 


M 


1 


4178 4 


3.89 x 10 3 i 


M 




1.5165 4 


5 x 10~ 4 i 


H 


1 


4178 4 


3.90 x 10 3 i 


H 




1.5188 4 


5 x 10" 4 i 


A 


1 


4196 4 


3.93 x 10- 3 i 


A 



Table 2: Experiments on the obstructed waveguide. (Accuracy: S, silly; L, low; M, medium; H, 
high; A, from Aslanyan et al. [AsPaVaj.) 
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The agreement with the results of Aslanyan et al. is good for the case of obstacle displace- 
ment 5 = (the case of an embedded eigenvalue). In the cases 5 > the agreement is less 
good, but the error is no worse than 2 x 10~ 3 . The values which we calculate change very little 
between the medium and high accuracy cases. 



15 Generalizations — resonances of a scatterer 



I n § 14 we described how to calculate resonances in the case when the domain has cylindrical 



ends. In fact a similar approach can be used to calculate resonances due to scattering by an 
obstacle. We consider a situation similar to that in Fig. [4] in which the domain is the region 
exterior to an obstacle (the shaded region). The cylindrical ends are now replaced by the region 
x 2 + y 2 > R 2 and the boundary T is the circle of centre (0, 0), radius R. 




Figure 4: Scattering by an obstacle. The transverse boundary T is replaced by a circle of radius 
R, the cylindrical ends by the region x 2 + y 2 > R 2 . 



As in § 14 we can compute resonances by replacing the map 1Z C X of (8.3) by the appropriate 
Neumann to Dirichlet map for the exterior region. A simple separation of variables shows that, 
in terms of appropriate Hankel functions H n , 



(15.1) 



in which the Fourier coefficients g n are 



2tt 



1 



g n = -= / «7(0)e-^#. 
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Suppose that we are interested in resonances which lie in the half-plane Im(A) < 0. Then 
choosing the branch of the square root so that Im(vA) < 0, the required Hankel functions are 
the ones which grow as 



(15.2) 



H n (r^X) 



^ r -l/2 e -rImVA 



as r — > +oo. These are the Hankel functions of the first kind. 



The matrix T of (11.5) is now replaced by a matrix which is clearly still diagonal. If we 
order the functions {e in& \ n £ Z} in the order 



■} 



then the corresponding matrix T is 
(15-3) 

_ .. / H (R^/X) h^rVx) H-i(RV\) h 2 (rVx) h_ 2 (rV\) 

T = diag 



OHS(WA) VXH[(RVX) VXH'^RVX) VXH 2 (RV\) VXH'_ 2 (RVX) 
The resonances can then be calculated following the same procedure as in § [14| 

15.1 Example 3: a cavity resonance problem 

We consider scattering by a C-shaped barrier almost enclosing a cavity as shown in Fig. [5] 




Figure 5: Cavity Resonator. The points P, Q, R and S have coordinates (1, e), (1.1, e), (1.1, — e) 
and (1,— e) respectively, and are joined by two circular arcs 



On the boundary of the domain we impose Dirichlet boundary conditions. For the Laplacian 
on this domain, a ess = [0, +oo). By simple separation of variables it is easy to show that no 
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e 


= 0.2 




e = 


= 0.3 


medium 


5.5158 - 
13.8574 


4.8 x 10" 
-4 x 10" 


- 4 i 
H 


5.1987- 
12.7505 


- 0.0049i 

- 0.1045i 


high 


5.5141 - 
13.8554- 


4.6 x 10" 
- 9.7 x 10 


- 4 i 
- 3 i 


5.1976- 
12.8165 


- 0.0051i 

- 0.1025i 



Table 3: Resonance calculations for C-shaped resonator. (Levels of accuracy: 'medium' means 
4 mesh refinements, R = 1.5 and A max = 50; 'high' means 5 mesh refinements, R = 1.3 and 

Amax = 100.) 

non-trivial solution of the Poisson equation —An = An on this exterior domain can lie in L2 
for A > 0, no matter what boundary conditions it may satisfy. Thus there can be no embedded 
eigenvalues. 

Nevertheless it is known (see [BroHiMa]) that for small e this problem possesses a resonance 
close to each Dirichlet eigenvalue for the Laplacian on the unit disc. We wish to compute these 
resonances. 

We first create an 'inner domain' by introducing an artificial outer boundary which is a circle 
centred at the origin with radius R > 1.1. It is important that R not be too large: as R — > 00, 
the (Neumann and Dirichlet) eigenvalues of the inner domain become dense in any compact 
subset of [0,oo), and so in particular it becomes very difficult to distinguish resonances close to 
the real axis from Neumann eigenvalues of the inner domain. For numerical reasons it is's also 
important that R not be too close to 1.1: otherwise the annular region 1.1 < a/x 2 + y 2 < R 
becomes very narrow, requiring a very fine mesh, and the numerical calculations require too much 
computing time. We carried out our calculations for R = 1.3 and, as a check, for R = 1.5. 

The Dirichlet problem on the unit disc has eigenvalues 5.783 (simple) and 14.682 (multiplicity 
2). We started by calculating the first channel resonances of the C-shape resonator closest to 
these eigenvalues. The results are shown in Tablefjj The resonance closest to the first of these 
eigenvalues is rather difficult to calculate accurately because it is so close to a pole of 1Z° X , the 
Neumann-to-Dirichlet map of the interior domain whose poles are the Neumann eigenvalues of 
the interior domain. See Fig. [6| which shows a contour plot of 

|det(L inner (A)Lo U 1 ter (A) - 1)1 , 



where, in view of (14.1), L mrier is given by 



£inner(A) = Rl : M,l:Af + $l:M,l:N {Dl:iV,l:Jv(A) — Di : jv,l:jv(Arj)} {Si-M,1:N)* , 

while L outer (A) is the analytic continuation of Ti:M,i-.m{^) already described. The point close 
to 5.2 — 5 x 10~ 3 i is the zero of the determinant, and is the approximate resonance. It is not 
easy to see this zero so close to the nearby pole. This requires sampling of the determinant on 
a very fine grid, which is expensive. 

The resonance close to the second eigenvalue is not so close as to cause the same numerical 
difficulties. 

We also sought to calculate resonances further out in the complex plane. These are rather 
more difficult to obtain — see Fig. [7J where we compare the results for the medium and high 
levels of accuracy. The agreement is rather poor. In fact the change in R is more responsible for 
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5.17 5.175 5.18 5.185 5.19 5.195 5.2 5.205 5.21 

Re(k) 



Figure 6: Resonance very close to a Neumann eigenvalue for the interior domain; R = 1.5, 
e = 0.3. At lower resolutions these are virtually indistinguishable. This is a problem, as not all 
Neumann eigenvalues are close to resonances! The contours are the level sets of a determinant 
which is zero at the resonance 5.199 — 5 x 10 _3 i. Unfortunately the same determinant also has 
a pole at the nearby Neumann eigenvalue 5.187. 



this than the level of mesh refinement or the value of A max for fixed R. This is rather different 
from the situation in [Ma], where the resonances were calculated by complex scaling. There, 
the sensitive resonances depended more on the mesh than on the domain of truncation. 

15.2 Example 4: Scattering by a Gaussian potential 

We consider the resonances of the Schrodinger operator in L2(M 2 ) given by — A + q(x,y) where 
q is a superposition of three Gaussians: 

3 

q{x,y) = C^exp(-i/(x - xj) 2 - v{y - yj) 2 ), 

with (xx,yi) = (0,-1), (>2,y2) = (sin(7r/3),cos(vr/3)), (x 3 ,y 3 ) = (-2:2,2/2), C = 40 and 
v = 2. For large C there is a high potential barrier on the unit circle x 2 + y 2 = 1 which gives 
rise to almost-trapped modes. However because q is rapidly decaying the PDE — Au + qu = Xu 
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Resonances far from real axis (medium accuracy) Resonances far from real axis (high accuracy) 




Re(X) Re(X) 

Figure 7: Resonances far from the real axis for e = 0.2: the medium and high accuracy results 
are only qualitatively similar 



number of mesh refines 


R 


resonance 


2 


4.0 


4.257 - 5 x 10" 4 i 


3 


4.0 


4.140-5 x 10~ 4 i 


3 


5.0 


4.185 — 10- 3 i 


4 


6.0 


4.402 - 4 x 10" 4 i 



Table 4: Resonance calculations for Gaussian potential. (A max = 50) 

in fact has no non-trivial L2 solutions for A > 0. The almost-trapped modes are the real parts 
of resonances close to the real axis, which one may attempt to compute by the approach in 
Example 3. 

An inner domain is created with boundary on the circle x 2 + y 2 = R 2 . The value of R 
should be chosen so that for x 2 + y 2 > R 2 , the function q is small enough to be neglected, 
which allows the PDE to be approximated by the Helmholtz equation. This allows us to use the 
same expression for the outer Neumann-to-Dirichlet map as for Example 3. 

As in Example 3, resonances far from the real axis are very unstable and cannot be calculated 
reliably. However we were able to calculate some resonances close to the real axis with moderate 
accuracy. One case is listed in Table [4| 

Lin [Lij attempts to find these resonances using the complex scaling method and a specially 
developed variational PDE-solver with a particularly chosen set of basis functions. Even then, 
many spurious results are generated, and Lin claims to be able to filter them out. 
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16 Conclusions 

Although our original motivation for doing this work lay in the need to devise a numerical method 
which would be easy to implement, the numerical results show that the accuracy achieved by 
our approach is actually quite competitive with conventional techniques. Moreover, although 
Neumann series converge rather slowly, two features ensure that for these problems the run- 



times were not excessive: firstly, the fact that the quantities in the eqn. (11.2) are calculated 



once at the outset, so the most expensive calculations need not be repeated for each different A 



at which the Neumann-to-Dirichlet maps are evaluated; secondly, the trick in eqn. (13.4) acts 
as a convergence acceleration technique for the Neumann expansion and means that in many 
cases we are able to get away with very few eigenfunctions. 

There are further modifications which could easily be made to improve efficiency. Foremost 
among these would be to replace contour-plotting of the condition number k(R(X) + T(X)) by 
a routine based on contour integration for finding zeros of an analytic function. There are good 
practical reasons, however, for not doing this. Resonances close to the real axis can already be 
located quite accurately by finding the maxima, along the real axis, of k(R(X) + T(X)); this is 
easy since it is a one-dimensional search. As already observed by Abramov, Aslanyan and Davies 
[AbAsDa| . resonances further from the real axis are typically very unstable; to spend a lot of 
time accurately locating the resonances of an approximating problem is therefore pointless. 
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